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Abstract 

Many problems in physics are inherently of multi-scale nature. The issues of MHD turbulence or magnetic recon- 
nection, namely in the hot and sparse, almost collision-less astrophysical plasmas, can stand as clear examples. The 
Finite Element Method (FEM) with adaptive gridding appears to be the appropriate numerical implementation for 
handling the broad range of scales contained in such high Lundquist-number MHD problems. In spite the FEM is 
now routinely used in engineering practice in solid-state and fluid dynamics, its usage for MHD simulations has re- 
cently only begun and only few implementations exist so far. In this paper we present our MHD solver based on the 
Least-Square FEM (LSFEM) formulation. We describe the transformation of the MHD equations into form required 
for finding the LSFEM functional and some practical issues in implementation of the method. The algorithm was 
tested on selected problems of ideal (non-resistive) and resistive MHD. The tests show the usability of LSFEM for 
solving MHD equations. 
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1. Introduction 

Dynamics of magnetized plasma at sufficiently larg e spatial and temporal scales can be adequately described by 



the set of magnetohydrodynamic (MHD) equations lllTll . In many problems we face the situation with high Lundquist 
(a.k.a. magnetic Reynolds) number 

S = Reu = LfioVA/ri , 

where L is the characteristic size of the system, Va - Bj yfp the typical Alfven velocity {B and p being the mag- 
netic field strength and plasma density, respectively) and 77 the electric resistivity. A direct consequence of the high 
Lundquist number is a large separation between the system size and the dissipation scale. The cascading fragmenta- 
tion of the current layer in the magnetic reconnection in solar flares can serve as an example of such a multi-scale 
problem: The span between the eruption size 10^ km) and the dissipation scale (1 m - 10 m) in the practically 
collision-less coronal plasmas easily extends seven orders of magnitude. 

In general, there are two approaches how to handle such a broad range of scales. The first one uses a moderate 
numerical resolution and models the physics on the sub-grid (unresolved) scales using some plausible assumptions on 
the micro-scale statistical properties (correlations) of the quantities that define the system (e.g. flow or magnetic field). 



Among them, e.g., the Large-Eddy Simulations (LES) 01411 or Reynolds-Averaged Numerical Simulations (RANS) 



01211 belong to the well known methods used widely in engineering applications in the fluid dynamics. 

The second approach is based on direct simulations that cover all the scales contained in the problem. Tradi- 
tionally, the Adaptive Mesh Refinement (AMR) technique is used with the Finite-Diff'erence/Finite Volume Methods 
in order to resolve high-gradient regions locally, keeping the total number of grid points required for simulation at 



a manageable level [iJ, |7|, |22|] . Nevertheless, also this approach has its limitations caused by introduction of arti- 



ficial boundaries between fine and coarse meshes. This problem, however, can be cured by the methods based on 
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unstructured mesh, such as is used in FEM. With this in mind we have implemented a FEM-based solver for MHD 
equations and present it in the current paper From various FEM formulations we have chosen the LSFEM because 
it is robust, universal (it can solve all kinds of partial differential equations) and it is efficient - it always leads to the 
system of linearized equations with symmetric, positive definite matrix OlOll . The LSFEM keeps many key properties 
of the Rayleigh-Ritz formulation even for systems of equations for which the equivalent optimization problem (in 
Rayleigh-Ritz sense) does not exist |i4J. 

Despite of the FEM applications in the fluid dynamics made a substantial development in the past years, its usage 
for numerical solution of MHD equations is still rather rare. For example, the NIMROD f^T] and M3D codes 190 - 
based on Galerkin formulation - belong to a few known implementations of FEM-based MHD solvers. Related work 
also has been done by Lukin 01 311 who implemented the MHD (and two-fluid) equations within the more general code 
framework SEL based on the Galerkin formulation with high-order Jacobi polynomials as the basis functions. 
However, to our knowledge, the LSFEM implementation of the MHD solver described in the current paper is the first 
attempt of this kind. 

The paper is organized as follows: First, we briefly describe the underlying MHD model. Then, the MHD equa- 
tions are re-formulated in the general flux/source (conservative) formulation. Temporal discretization, reduction to 
the first-order system, and linearization procedure are described subsequently. Then, the properties of the least-square 
formulation of FEM are briefly summarized. Some practical arrangements of the LSFEM implementation of the MHD 
solver follows. Finally, the code is tested on a couple of standardized model problems and the results are discussed 
with respect to the intended application of the code to the current-layer filamentation and decay during the magnetic 
reconnection in solar eruptions. 



2. MHD equations 

The large-scale dynamics of magnetized plasma can be described by MHD equations for compressible resistive 
fluid d: 

d,p + V-(pM) = 

pdtU + p(u • V)h - -Vp +j X B + pg (1) 
d,B = Vx(h X B) - Vx(i]j) 
d,U + V -S =pu- g, 

where p, v, B, U are density, macroscopic velocity, magnetic field, and total energy density, respectively, g being the 
gravity acceleration. The energy flux S and auxiliary variables j (current density) and p (plasma pressure) are given 
by the following relations: 

V X B = poj 

U^-^ + -pu^ + — (2) 
y-l 2 2/io 

U + P+ M - -^^ -B + —jxB 

2po j Pa Po 

In the (almost) collision-less plasma, in which we are mostly interested, the classical resistivity usually plays a 
small role. Instead of that various microscopical (kinetic) effects influence the plasma dynamics via other terms in the 
generalized Ohms law |5]. In order to mimic these processes, whose modeling is beyond the scope of MHD approach, 
we re-consider the parameter 77 as a generalized resistivity, including the eff'ects like wave-particle interactions or off- 
diagonal components in the electron pressure tensor into it. As such effects are - in general - observed in the highly 
filamented, intense current sheets we model the anomalous generalized resistivity as follows: 

{0 : Iv'dI < 

L . |V'dI > Vcr 
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Thus, the non-ideal effects are turned on whenever the current-carrier drift velocity 



VD(r, t) 



\jir,t)\ 



(4) 



exceeds the critical threshold Vcr- 

In order to solve the Eqs. ([1) numerically, it is convenient to rescale all the quantities into the dimensionless units. 
Thus, all the spatial coordinates are expressed in the characteristic size L and times in Alfven transit time ta - L/Va, 
where Va = Bq/ sjpo is a typical Alfven speed. Magnetic field strength B and plasma densit y p are given in units 
of their characteristic values Bo and po and similar scaling holds for the other quantities - see 111 ill or B for details. 
From now on we shall use this dimension-less system. 

In order to utilize a more universal LSFEM implementation for more general form of equations [c.f. with SEL 
approachT?) the set of MHD Eqs. ([U is rewritten into the conservative (flux/source) formulation: 



d,'¥ + 5,- f /CP, d,.'¥) = S(xj, «P, f) 



(5) 



Here the local state vector ^ - (p,n, B, U), n - pv being the momentum density. The flux F and the source-term S 
are defined as 



(6) 



n 




f ^ 1 


pvv -BB+ iixiip + B^) 


, s = 


Pg 


^3x3 • E 







, {h + Ek)v + 2ExB 







where /sxs is the 3x3 unit matrix, 63x3 is the permutation pseudo-tensor, E - -v X B + r]j is the electric field strength. 
The the enthalpy and kinetic energy densities are h - ypliy - 1) and Ek - pv^, respectively. 

3. FEM formulation of MHD system 

In general, FEM is formulated for the linear problem 

Lh = / inQ 



(7) 



Bh = g onF = c?Q 



where L is the linear (diff^erential) operator, B the boundary operator, Q is the domain and F is the boundary of Q. 

In order to reformulate the system of Eqs. (|5) into the standardized problem (|7} several steps have to be undertaken. 
First of all, we perform the time discretization. We use the standard 0-differencing scheme [see, e.g.,01: 



(8) 



where parameter © e (0, 1) controls the implicitness of the scheme, and n and n + 1 designate the old and new 
time-steps, respectively. The scheme leads to the following semi-implicit equation 



(9) 



where the RHS vector R" consists of components known at old time step. 

Since practical implementations of LSFEM require first order system of PDEs fl we further transform the 
system (|5]l to the required form introducing a new independent system variable - the electric field 



-V X B + 77V X B 



(10) 



The procedure is basically analogous to the velocity-vorticity formulation of the Navier-Stokes equations in the CFD. 

3 



A frequent problem in the numerical MHD is a violence of the solenoidal condition pn, = V • B = 0, where 
the (dummy) variable represents the artificial density of the magnetic charge. The advantage of the LSFEM 
implementation is that this constraint can be directly included into the set of the governing equations [10]. Then 
assembling the solenoidal condition together with Eqs. (|9]i and ( fTOl i we arrive to the following 1 st-order vector equation 
for our modeled system: 





d 

dxi 


' TFi('¥,E) ' 




' tS(T) ^ 
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(11) 



where all the LHS terms are evaluated in the advanced time-step « + 1 . Here, F(*P, E) and 5(*P) are given by Eq. (|6]l 
with E considered as an independent variable now and t = &At. The fluxes G = — f3x3 • B and H - -B imply from 
Eq. ([Tol l and the solenoidal condition. The source term component r(*P) = {nB - Bn)lp - 6(4*) • V;/. We keep the 
artificial magnetic-charge density p,n at zero. 

Eq. (fTTT l can be written in the conservative form similarly as in Eq. (|5]l 



* + 5,. F,^*) - 5(*) = R 
with the extended state vector, fluxes and source terms in the form 



(12) 
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(13) 



Since the extended flux F and source term S depend non-linearly on the state vector a linearization procedure 
has to be applied in order to transform the svste m ([TZ j into the FEM-conforming form ([7). We use the standard 
Newton-Raphson (NR) iterations in each time step lllOi 12011 . Thus, the flux at the NR iteration k + 1 can be expressed 
as 10]: 

f ,.(T*^i) = Fii¥) + 4 -("f'^' - ¥) (14) 
and analogous expression holds for the source term. Introducing the Jacobians 



A, = d^Fi 



c = d^s\^ 



(15) 



the final equation for NR iterations reads 



R + (5,.,A/ + C) •¥ - S(¥) , 



(16) 



where the RHS contains only the terms from the k-th iteration of the currently solved time-step n + 1 and variables 
known at the previous step n. Eq. (fTSI l is already in the form (|7) with 



L = \t + d,Ai + Ai^+C 

OX; 



f^R + {d,Ai + C) ¥ - S(¥) 



(17) 



4. LSFEM implementation 

In the least-square formulation of the FEM the problem described by Eqs. Q is transformed to seeks the minimum 
of the functional 

(m) = j (Lu- ff AQ. + W J (Bm 



gfdr 



(18) 
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where w is appropriate mesh-dependent weighting factor ||4|]. As in other FEM implementations, the solution is 
searched for in a limited subspace of functions that are formed as a union of the piece-wise functions Ug{x) defined on 
a single, in our code triangular element, as a linear combination of the basis functions <i>i{xy. 

u,(x) = ul%(x) , (19) 

where u'^ can always represent the value of a function u(x'j) in a properly selected point (the node) x'j. Here / denotes 
element-wise index of the node. In our code we use Lagrangian polynomials for basis functions ^/(x). 

Varying the functional (fTSl l and inserting the expansion (fT9] ) we arrive to a set of linear algebraic equations for 
each internal element in the form 



J]i^J L^<l),L<DjdQj ■ = J L^<[),/dQ, (20) 



where fig c f2 is the domain of the e-th element in the global domain Q. The boundary elements contain additional 
terms obtained from the boundary operator (the second term in Eq. (fT 8] l). For fast evaluation of local integrals we use 
Gaussian quadrature |0] in the system of element natural coordinates jlSll . 

Equations (l20l i for each element are finally assembled to a global linear system of equations via mapping the 
element- wise node index j to a global node index described in 16J, to obtain 

Y^KiN-u^^f. (21) 

N 

The final matrix K is sparse, symmetric and positive definite. In our code we use preconditioned Jacobi Conjugate 
Gradient Method (JCGM) |T6^ for solution of the system (1211 1. 
The entire algorithm can be summarized as follows: 

- time loop - adapt time step size according to CFL condition, check final desired time 

- Unearization loop - if ' ^^t^,^ ' < e or maximum iteration count is reached continue to next time step 

- assembling stiffness matrix K element by element 

- integration by Gaussian quadrature 

1 . compute the operator matrices for each basis function 

2. multiply the operator matrices then add the result into stiffness matrix 

3. multiply the operator matrix by the RHS then add result into the load vector 

- next Gaussian point 

- next element 

- find new solution m*^^' of system (ISTT i by the JCGM 

- next linearization 

- next time step 

Thanks to the iterative nature of the JCGM, the solver algorithm can be rather easily parallelized via MPI. We 
decompose the entire domain into subdomains, splitting the global matrix K and the load vector / into corresponding 
segments with rather small overlaps related to internal-boundary nodes shared by both adjacent subdomains. Matrix 
multiplications are then performed only locally (per-process) and, finally, resulting global vectors are appropriately 
assembled using MPI operations that transfer the data related to overlapping nodes only. 

5. Numerical Tests 

In order to assess usability and properties of the LSFEM MHD solver we perform several tests on standardized 
ideal (non-resistive) and resistive MHD problems. For all test we use the adiabatic index y - 5/3, the implicitness 
parameter = 0.5 (Crank-Nicholson time discretization), and the Courant number 0.6. 




Figure 1: The LSFEM solution of the MHD shock tube test (Ryu- Jones problem) at time f = 0. 1 with the first-order basis functions, (a) density 
profile (red dashed line) and By profile along A-axis. (b) t', profile along ,t-axis. 




Figure 2: The LSFEM solution of the MHD shock tube test (Ryu- Jones problem) at time f = 0.1 with the second-order basis functions. Displayed 
quantities are the same as in Fig.^ 



5.7. Ryu-Jones discontinuity test problem 

First, we applied our code onto the standard Ryu-Jones ideal MHD ID shock/discontinuity problem lfl9ll . The ini- 
tial state is given by prescriptions (p, Vi, v,,, v~, B^, B,, B,, £) = (1, -1, 0, 0, 0, 1, 5, l)in the left half, and (p, y^, v,,, v^, B^, By, B,, E) = 
(1, 1, 0, 0, 0, 1, 5, 0) in the right half of the computational box, respectively. The domain (-0.5, 0.5) was divided into 
512 elements. We used the first order basis functions to approximate the FEM solution. The boundary conditions 
on both ends are of von Neumann type. Results at time f = 0.1 are shown in Fig. [1] They correspond and could be 
compared with Fig. 3b in ifioll . 

In order to study influence of basis-function order on the approximate solution we calculate the same test problem, 
now with the second-order Lagrange polynomials. All other parameters are the same as in the previous case displayed 
in Fig.[T] The results are shown in Fig.|2] 

5.2. Orszag-Tang vortex test problem 

A next test we performed standard Orszag-Tang 2D ideal-MHD vortex problem l,15J . The initial state was given 
by the following relations: 



p ^ py 
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Figure 3: The Orzsag-Tang vortex. The color coded plasma density is displayed at time t = 0.25 (a) and t = 0.50 (b). 



= 

P = 



• sin(27r3;), 
1 



TTy — sin(27rji:) 



V47r 



sin(27r3'), B,, = 



1 



V4^ 



sin(47rjic) 



The computational domain 1 .0 x 1 .0 was discretized by 2 x 640 x 640 triangular elements. We apply periodic boundary 
conditions at all boundaries. The first-order basis functions were used in this simulation. Results in Figs. [3]and|4]show 
the plasma density and the magnitude of the magnetic field, respectively, at times t - 0.25 (a), and t - 0.50 (b). 

5.3. Resistive decay of a cylindric current 

In order to assess the applicability of our code to the solutions of non-ideal (resistive) MHD problems and to 
estimate its numerical resistivity we performed a following test: At the initial state f = a cylindrical current j{r) - 
(0,0,7,(r))with 

;oJo(jca'^) : r < r^ 
■■ r> rf, 



jz(r) 



(22) 



flows through a plasma of a uniform density po. Here, jo = 1 is the amplitude of the current density on the cylinder 
axis, ro = 1 is the cylinder radius, and ~ 2.40 is the first null of the Bessel function of the 0th order Jo(x). The 
resistivity inside the cylinder (r < ro) is uniform 77 = 770 = 0.1, outside t] - Q. In order to be able to compare 
the numerical results with an analytical solution and to split advective and resistive properties of the code we set all 
velocities to zero at f = and the density to a very high value po = 10^ to keep the plasma in rest. In the limit p — » 00 
the MHD system ([T]) effectively reduces into the diffusion equation 



d,B + = 

whose analytical solution for our initial state keeps the form j(r, t) - (0, 0, t)) with 

jzi'; t) = johixN—) exp(-7f) + -^Ji (xa,) [1 - exp(-yf)] S{r - ro) 



(23) 



(24) 




Figure 4: The Orzsag-Tang vortex. The color coded magnitude of the magnetic field is displayed at times t = 0.25 (a) and ( = 0.50 (b). 

where Ji(ji:) is the Bessel function of the 1st order and 5{x) is the Dirac deha function. The decrement y reads 

The second term in Eq. (l24l l represents an induced surface current that compensates resistive decrease of the current 
density inside the column to keep the magnetic field in the outer super-conducting domain constant. The corresponding 
magnetic field is of the form B - (0, B^, 0) where 



for internal (r < ro) region and 



ro r 
B^{r,f) = ;■() — Ji(x/v— )exp(-yf) 

Xn ro 



B^(r,t) = jo—h(xN) 

Xn 



(26) 



(27) 



for the outer space. 

Computational domain is divided into a homogeneous mesh of 2x512x512 triangles in our numerical test. We 
use the first order basis functions to approximate the numerical solution. Free boundary conditions were applied on 
all boundaries. The results of this test are shown in Fig. |5] Fig. |3a) shows time evolution of the current density profile 
along y = for five subsequent time instants. Resistive decrease of j. inside the column accompanied by formation 
of the induced surface current are well visible. Fig. |5}b) shows a comparison of numerical and analytical solutions for 
time evolution of the current density jz(x,y, t) at x - Q,y - Q. 



6. Discussion and conclusions 

The FEM represent an alternative to FDM/FVM that are traditionally used for solution of MHD problems in 
astrophysics. Its attractivity implies from its unstructured mesh that allows for appropriate local refinement without 
formation of qualitative internal boundaries between the fine and coarse meshes. This property makes it very useful 
for handling the multi-scale problems, for example the problem of magnetic reconnection in solar flares ||2| (and other 
large-scale systems) or MHD turbulence. 
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Figure 5: A resistive decay of a cylindrical cun'ent density with time, (a) profiles of j-{x, 0, t) at five subsequent times, (b) The time profile of 
j.(0,0,t) - comparison of numerical and analytical [Eq. | |24H solutions. 



With this intention in mind we have developed the LSFEM implementation of a MHD solver whose descriptions 
and preliminary results from its application to the standardized test problems are presented in this paper 

To sum up the main points of our implementation: Transformation of the MHD equations ^ to the standard 
FEM problem (Q involves several steps: (i) Standard 0- time discretization, (ii) Decrease of the order of the system 
of equations by introduction of a new variable - electric field strength, and (iii) Newton-Raphson linearization. The 
possibility to include the solenoidal condition V • B = directly into the system of equations certainly belongs to 
advantages of LSFEM formulation, as well as a natural involvement of the boundary conditions. The element-by- 
element assembling of the global stiffness matrix and the iterative nature of JCGM solver allow for rather easy and 
efficient MPI parallelization. Integrals over elements are efficiently performed via Gauss quadrature. 

We performed several standardized tests focused on an ideal and resistive MHD. The LSFEM MHD solver quite 
closely reproduces results published for the Ryu-Jones shock tube problem 1 19]. Small spurious oscillations appear 
around the points where the first derivative of an analytical solution does not exist. Choice of the higher-order basis 
functions makes the situation even slightly worse. 

Similar feature can be seen in the results from the Orszag-Tang vortex test problem. While the large-scale dy- 
namics agree well with those obtained from the 'gauge' codes, small oscillations accompanying the shocks are visible 
again. These effects are caused by the least squares curve fitting approach [4]. We believe that it can be cured by 
an introduction of a small background resistivity and local refinement of the mesh around the discontinuities, with 
the element size corresponding to the resistivity-controlled (magnetic) Reynolds number Such approach is fully in 
line with the intended usage of the code for detailed studies of the current sheet filamentation and fragmentation in a 
large-scale magnetic reconnection in solar flares. Indeed, in the solar corona we have a very small background classi- 
cal resistivity due to (rare) collisions between electrons and ions as well. Hence, having mesh around the filamenting 
current sheet locally refined as much as possible we can set the background (physical) resistivity accordingly and 
approach thus the realistic Lundquist number in the solar corona. 

Finally - with the intended usage of the code in mind - we have tested the properties of our implementation 
for solution of the resistive problems. In order to get a comparison with an analytical solution we have 'frozen' the 
plasma dynamics by setting high matter density and we concentrated on a purely diffusive problem. The results show 
a rather good agreement with the analytical solution. Namely, the induced surface current density is located only at 
a few elements and did not diffuse further with time. This is an important result for intended future studies of the 
current-sheet filamentation in the flare reconnection. 

The tests show basic applicability of our LSFEM implementation of the MHD solver for a solution of selected 
problems. At the same moment they reveal the necessity to involve both the adaptive spatial refinement (it has already 
been implemented) and adaptive change of the order of basis functions over selected elements (h-p refinement). These 
features will be implemented into our code in a near future. 
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